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We study the thermodynamics of the two-dimensional Hubbard model within the dynamical 
cluster approximation. We use continuous time quantum Monte Carlo as a cluster solver to avoid 
the systematic error which complicates the calculation of the entropy and potential energy (double 
occupancy). We find that at a critical filling, there is a pronounced peak in the entropy divided 
by temperature, S/T, and in the normalized double occupancy as a function of doping. At this 
filling, we find that specific heat divided by temperature, C/T, increases strongly with decreasing 
temperature and kinetic and potential energies vary like In T. These are all characteristics of 
quantum critical behavior. 



Introduction- The properties of the hole-doped 
cuprate phase diagram, including a pseudogap (PG) at 
low doping and unusual metallic behavior at higher dop- 
ing, have led many investigators to propose that there is 
a quantum critical point (QCP) in the cuprate phase di- 
agram at optimal dopingi. Some of the most compelling 
evidence for the QCP is from various thermodynamic ex- 
periments^. 

In the previous work^ employing the dynamical cluster 
quantum Monte Carlo (QMC) method^"^ to calculate the 
quasiparticle fraction for the two-dimensional (2D) Hub- 
bard model, we found evidence for a QCP separating a 
region of Fermi liquid (PL) character at high doping from 
a region of non-Fermi liquid PG character at low doping. 
At the critical doping, we found marginal Fermi liquid 
(MFL) character— which is also seen above the PL and 
PG temperatures. 

In this paper, we provide further evidence for a QCP 
in the thermodynamic properties, including the ener- 
gies, the specific heat and entropy. At the critical fill- 
ing, both the kinetic and potential energy data show a 
In T low-temperature dependence, leading to a T In T 
specific heat. The entropy divided by temperature, S/T, 
is strongly peaked at the critical point. In contrast to 
the results of the previous studj*^, these quantities are 
independent of the location and character of the Fermi 
surface and thus avoid ambiguity in their interpretation. 

Formalism- We start with the 2D Hubbard Hamilto- 



nian 



(1) 



where c|^^(ckcr) is the creation (annihilation) operator for 



electrons of wavevector k and spin a 



,c,vr is the 



number operator, = — 2t (cos(fc2;) -|- cos(fcy) — 2) with 
t being the hopping amplitude between nearest neighbor 
sites, and U is the on-site Coulomb repulsion. The pri- 
mary focus of this Rapid Communication is the behavior 
of the energies. For lack of a better term, we will refer 
to the quadratic part of Hk as the kinetic energy and the 
interaction part as the potential energy. The correspond- 



ing energies per site are^ 



E, 



N N 



(2) 



^ ^-^^^ S.(k,zc.„)a(k,zc.„), (3) 



N 2N 

where G'(j(k, ia;„) and Scr(k, are the single-particle 
Green function and self-energy respectively, with w„ = 
(2?7 + 1)itT and represents the number of sites. 

Aside from the numerical calculations, we obtained an- 
alytically the leading low-T behavior for Ek and Ep for 
each of the regions around the QCP starting from the cor- 
responding forms for the self-energy. In the FL region, 
the imaginary part of the self-energy has the form 



-apL max (lj^, T^) , \uj\ < ujx, 



(4) 



where lox is a cutoff frequency and apL > 0. In the MFL 
region, this quantity becomes^ 



^'mfl (^) 



-aM_FLniax(|w| ,T) , |a;| < (5) 



where olmfl > 0, and tOc is a frequency cutoff. Finally 
in the PG region, we consider the following Ansatz- for 
the imaginary part of the self-energy: 



(6) 



where ui* is a cutoff frequency and A (k) = 
Ao (cosfc^ — cosfcj,) with the constant Aq being the PG 
magnitude. 

Inserting the forms (gl), ^ and © in ^ and © and 
performing a low-T expansion in the resulting integrals, 
we obtained that both Ep and Ek exhibit a leading 
low temperature behavior in the PG and FL regions. 
Both energies display a In T behavior consistent with 
a In r total energy^ in the MFL region. 

Methodology- We solve the Hubbard model within the 
dynamical cluster approximation (DCA)^. The DCA is a 
cluster mean-field theory which maps the original lattice 
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model onto a periodic cluster of size Nc = embed- 
ded in a self-consistent host. Spatial correlations up to 
a range Lc are treated explicitly, while those at longer 
length scales are described at the mean-field level. How- 
ever the correlations in time, essential for quantum criti- 
cality, are treated explicitly for all cluster sizes. To solve 
the cluster problem we use weak coupling expansion con- 
tinuous time QMC (CTQMC) methodi^iii with highly 
optimized blocked and delayed updates^^. In the previ- 
ous work^, a Hirsch-Fye QMC (HFQMC) algorithm was 
used. However, the Trotter error intrinsic to HFQMC, 
which is particularly large for the potential energy (dou- 
ble occupancy), prevented us from calculating S since 
systematic errors from multiple temperatures tend to ac- 
cumulatei^. The CTQMC algorithm eliminates the sys- 
tematic error, and generally has a lower overall compu- 
tational cost than HFQMC, making it a better choice for 
these calculations. Unless otherwise displayed, the sta- 
tistical error bars in these calculations are smaller than 
the plotting symbols or line widths used. 

Here, we study the normal state on a A^c = 4 x 4 site 
cluster for U = 6t. We confirmed that increasing the 
cluster size to Nc = 24 sites does not significantly alter 
our results. We obtain the energies from Eqs. ^ and 
dH]) , and, following Ref. [H, the entropy from 

S{f3,n) = S{0,n)+PE{/3,n)- E{f3',n)dp', (7) 



where S'(0, n) = -n In f - (2 - n) In (l - f ) , n is the fill- 
ing and f3 ~ 1/T. Finally, since the DCA preserves ther- 
modynamic consistency^^, our results obey the Maxwell 
relation 



dS_ 

dn 



T.1J 



9/1 
df 



(8) 



U,n 



where /i is the chemical potential. 

Results- Fig. [T] shows the normalized double occu- 
pancy, AD/ri^, as a function of temperature at fillings 
n = 1.00, 0.95, 0.85 and 0.70 together with fits to the 
MFL form for n = 0.95, 0.85 and 0.70, and the FL form 
for n = 0.70. For n = 0.95 and in a low temperature 
window, the Ep first increases with decreasing temper- 
ature and then reaches a maximum at a temperature 
which coincides with the PC temperature (T*) found in 
a previous study^. Further decreasing the temperature 
into the PC phase, Ep deviates from the MFL charac- 
ter and begins to decrease. This indicates that close to 
half- filling, where the interactions are more relevant, the 
system tries to gain energy by lowering the Ep. A more 
pronounced decrease in the Ep can be seen below the PG 
temperature in the half-filled case. At the critical filling, 
n = 0.85, Ep fits the MFL analytic form very well down 
to the lowest temperature reached. At this temperature, 
we find that the normalized double occupancy is maxi- 
mal at the critical filling (see the inset of Fig. [T]). For 
n = 0.70, the Ep increases with decreasing temperature 
as in the MFL region. However, we find that below the 
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FIG. 1: (color online) Temperature dependence of normal- 
ized double occupancy 4D/n^ . For a range of temperatures 
and values of filling, D{T) fits well the MFL form (dotted 
lines): D(T) = Do + aT^ ln(r/tJc) (Do, a and luc are fitting 
parameters). At low temperatures, D{T) shows FL behavior 
(dashed line) for n = 0.70 {D{T) = Do + aT'^, q < 0) and PG 
behavior for n > 0.85. Indicated by arrows is the crossover 
temperature T* — O.lSt for n = 0.70 along with the crossover 
temperatures Tfl = 0.24t and Tfl = 0.3t for n — 0.95 and 
n = 1.00 respectively. For n = 0.85, marginal Fermi liquid 
behavior persists down to lowest temperatures accessible with 
QMC. Inset: Normalized double occupancy versus filling for 
T = 0.077t, showing a peak at the critical filling n = 0.85. 



FL characteristic temperature (Tfl) (see Ref. Q), it de- 
viates from the MFL form and fits better to a quadratic 
function in T. Note that at intermediate temperatures 
(T* or Tfl <T < OM), the MFL form fits the Ep very 
well in all three doping regions. 

In Fig. [21 we show Ek versus temperature at the same 
fillings. Tfl and T* in Fig. [5] represent the onset of FL 
and PG regions obtained from Ek fits, and they have 
the same value as those shown in Fig. [T]for the Ep. At 
low temperatures, while the Ek deviates from its MFL 
fit for n = 1.00, 0.95 and 0.70, one sees no sign of a 
characteristic temperature scale at the critical filling. 

We note that the low-T behavior of both Ep and E^ fol- 
lows the analytical forms in the FL and MFL regions but 
slightly deviates from the predicted form in the PG 
region at low T. Also the energy scales below which the 
data deviate from the MFL T^ log T behavior vanishes 
at a critical filling, Uc = 0.85, where it persists down 
to the lowest temperature accessible, consistent with the 
existence of a QCP between FL and PG regions. 

The kinetic and potential energy data, together with 
the analytic forms of Eqs. ^ and ^ suggest that the 
total energy away from the half filling may be fit to the 
form 



E{T) = E{0) + Af{T)T^ + B{1~ f{T)) In 



(9) 



where /(T) = 1 / (cxp {{T - Tx) / 9) + 1) describes the 
crossover from the MFL to the quadratic behavior, char- 
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FIG. 2; (color online) Temperature dependence of kinetic en- 
ergy. MFL form fits data for a range of intermediate temper- 
atures (dashed lines). For n = 0.70, crossover to FL behavior 
is indicated as Tfl, while for n — 0.95 and n = 1.00, crossover 
to PG region at T* is evident by strong departure of kinetic 
energy from MFL fit, which for n = 1.00 even results in no- 
ticeable increase of kinetic energy with lowering temperature. 
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FIG. 3: (color online) Total energy per site E versus tem- 
perature for different fillings. The data are fit to a crossover 
form of the energy, Eq. ^ (dashed lines). The values of Tx 
determined from the fit are indicated as Tfl for n = 0.70 
and T* for n = 0.95. In the inset, the specific heat calculated 
from the fit is plotted versus temperature. 



acteristic of a FL or presumably PG region. A, B, 6, 
Tx and are the fitting parameters of the QMC energy 
data, as shown in Fig. [3l At low T, the fit is indistin- 
guishable from the data for all fiUings. The values of Tx 
determined from the fit are indicated as TpL ~ 0.1 5i for 
n = 0.70 and T* « 0.2At for n = 0.95 in agreement with 
the Ep and Ek fits. The calculation of the specific heat is 
a notoriously difficult problem and usually involves a fit 
of E{T) to a regularized (smooth) functional for m^^i^^ . 
Here, we already have an excellent fit, so wc obtain C/T 
from a derivative of the fit divided by temperature. The 
result is shown in the inset. As expected for n — 0.70, 
at low temperatures, C/T is fiat in T consistent with FL 
formation. This behavior is also seen for the data in the 
PG region, n = 0.95, but the data at the critical filling 
n = 0.85, shows a weak divergence at low T consistent 
with quantum critical behavior. 

The entropy per site exhibits intriguing behavior near 
the critical filling as the system is cooled, confirming 
much of the physics seen in C/T without the need for a fit 
or a numerical derivative. It is more strongly quenched 
in the FL and PG regions than in the MFL region with 
decreasing T, resulting in a maximum in S/T versus fill- 
ing at ri = 0.85 at the lowest accessible temperature (see 
Fig. [4] (a)). The fact that S/T continues to rise at the 
critical doping as T — > is consistent with the increase in 
C /T. Similarly, the low temperature S/T curves at dif- 
ferent temperatures nearly overlap for n < 0.85 in agree- 
ment with a constant C/T and indicative of a FL. One 
can see the strong decrease in S/T with decreasing tem- 
perature when the system crosses over to the PG phase 
around T* at large fillings (e.g. at half-filling). In fact in 
the PG region, n > 0.85, S/T coincide for different tem- 
peratures only at one filling, roughly n — 0.95. The flat 



low-T C/T seen at this flUing is accidental, and for higher 
filling one would expect C/T to be strongly suppressed 
at low T^. 

According to Eq. ([H]), a local maximum in S/T versus 
n corresponds to a flat chemical potential as a function of 
temperature. Therefore, at low T, the critical flUing can 
be identiflcd from the temperature dependence of /i for 
different fillings. This is shown in Fig. 3] (b) where, one 
can see that the near temperature independence of fi at 
n = 0.90 for 0.25t < T < 0.50t evolves into a broad max- 
imum centered around T = 0.15t for n = 0.87 which pre- 
sumably moves to n = 0.85 at low enough temperatures. 
These observations are consistent with the evolution of 
the maximum in S/T versus n as the temperature is low- 
ered from 0.50t to 0.08t (see Fig.|3](a)). In analogy to the 
half-filled case, a stationary chemical potential can be the 
signature of local particle-hole symmetry. This is consis- 
tent with the observation of near particle-hole symmetry 
in the cuprates in the proximity of optimal dopingi^. It 
is also in agreement with previous results showing that 
at n w 0.85 and for the same model parameters, the low 
energy density of states displays particle-hole symmetrj*^. 

Discussion- The results presented here strongly fa- 
vor an interpretation involving a QCP as opposed to a 
simple crossover from a FL to a non-Fermi liquid as the 
filling increases towards one. The In T behavior of the 
kinetic and potential energies, the peak in S/T which 
sharpens as T falls, and the logarithmic behavior of the 
speciflc heat C/T ~ InT at n = 0.85 together with the 
wide range of temperatures and fillings influenced by the 
critical point arc clear signatures of quantum criticality. 

The peak in the low temperature normalized double 
occupancy observed at the critical filling is also inter- 
esting. Near half filling AD/ii? is suppressed by strong 
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FIG. 4: (color online) Left panel: S/T filling dependence 
showing emergence of a peak at n = 0.85 at low temperatures. 
Right panel: Chemical potential temperature dependence for 
a range of fillings with PG and FL energy scales shown as 
T* and Tfl for n — 0.95 and n = 0.70, respectively. Note 
that the position of the maximum of entropy in the left panel 
corresponds to dfi/dT = in the right panel. As temperature 
is lowered, the maximum of entropy shifts towards lower n, 
causing local particle-hole symmetry for n — 0.85 at low T 
(see text). 



correlations. As the temperature is lowered {T < T*), 
it is suppressed further to allow the system to gain the 
most magnetic exchange energy. As the filling decreases 
towards zero, the Fermi energy will become smaller than 
U, and again AD/n^ is suppressed. So for any finite U 
we expect a peak in AD/n? at some finite doping. 

At the critical filling, in addition to the enhanced C/T, 



we find a peak in the charge susceptibility (not shown). 
These charge fiuctuations may be responsible for the lo- 
cation of the peak AD/v? near the critical doping. The 
peak in the charge susceptibility becomes sharper as the 
temperature is lowered, or when we include a next-near 
neighbor hopping t' > 0. This behavior suggests the pos- 
sibility that the QCP may be associated with a charge in- 
stability seen previouslyi^. This topic, as well as relation 
of the QPC and charge fluctuations to superconductivity, 
will be explored in an upcoming publication. 

Conclusion- We study the thermodynamics of the 2D 
Hubbard model using the dynamical cluster approxima- 
tion with CTQMC as a cluster solver. The latter elim- 
inates the Trotter error which complicates the calcula- 
tion of the entropy and potential energy. At the critical 
doping we find that C/T, obtained from a fit of E{T), 
increases strongly with decreasing temperature and the 
kinetic and potential energies are consistent with InT. 
Near the critical filling, we find a pronounced peak in S/T 
which grows as the temperature falls, consistent with the 
growth of C/T. We also find a peak in the normalized 
double occupancy as a function of doping. These are all 
characteristics of quantum criticality. 
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